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ABSTRACT 

The rates at which mass accumulates into protostellar cores can now be predicted in 
numerical simulations. Our purpose here is to develop methods to compare the statistical 
properties of the predicted protostars with the observable parameters. This requires (1) an 
evolutionary scheme to convert numerically-derived mass accretion rates into evolutionary 
tracks and (2) a technique to compare the tracks to the observed statistics of protostars. Here, 
we use a 3D-Kolmogorov-Smirnov test to quantitatively compare model evolutionary tracks 
and observations of Class protostars. 

We find that the wide range of accretion functions and timescales associated with gravo- 
turbulent simulations naturally overcome difficulties associated with schemes that use a fixed 
accretion pattern. This implies that the location of a protostar on an evolutionary track does 
not precisely determine the present age or final accrued mass. Rather, we find that predictions 
of the final mass for protostars from observed Tb i-Lb i values are uncertain by a factor of 
two and that the bolometric temperature is not always a reliable measure of the evolutionary 
stage. Furthermore, we constrain several parameters of the evolutionary scheme and estimate 
a lifetime of Class sources of 2-6- 10 4 yrs, which is related to the local free-fall time and 
thus to the local density at the onset of the collapse. Models with Mach numbers smaller than 
six are found to best explain the observational data. Generally, only a probability of 70 % was 
found that our models explain the current observations. This is caused by not well understood 
selection effects in the observational sample and the simplified assumptions in the models. 

Key words: Accretion, accretion discs - Methods: numerical - Methods: statistical - Stars: 
formation - Stars: statistics 



1 INTRODUCTION 

A major challenge still remaining in astronomy is to determine how 
mass accumulates to form stars. We now know from observations 
that almost the entire mass of a star li ke the Sun is gained during a 
short protostellar stage (Adams et al. 119871) ). Some of these proto- 
stellar objects are luminous yet cool, suggesting the existence of an 
even briefer phase within which the p rotosta r accretes rapidly from 
a highly obscuring core (Andre et al. Il993f0 . The existence of this 
so-called Class phase, however, highlights the need for the devel- 
opment of models which predict how a protostar abruptly evolves. 
A successful model would also predict the total accumulated mass 
or 'final mass', providing a fundamental input for the subsequent 
stellar evolution. However, it is debatable whether the development 
of low mass stars can be described by an evolutionary model or if 
the systems are simply too complex and diverse. 

To further the debate, we require three factors: reliable obser- 
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vations of protostars, models for how cores form and collapse out of 
a molecular cloud and a means of converting model parameters into 
observational parameters. Here, we develop a statistical method for 
comparing observational data to model data and proceed to test the 
method by making a specific comparison. This requires us to de- 
velop an appropriate technique to compare data within a three di- 
mensional parameter space. We first introduce, in turn, the data, the 
model and the interfacing scheme. 

The protostellar stage begins soon after the start of collapse 
of the cloud core. Half of the final stellar mass is gained in the 
very earliest, deeply embedded phase. Further accretion occurs dur- 
ing the Class 1 phase while an accretion disk still feeds material 
in the extended Class 2 stage. The number of observable Class 
sources at any given instant is relatively small, due to their brief 
lifespans. In addition they are difficult to detect since their mas- 
sive envelopes cause their spectral energy distributions (SED) to 
peak in the far infrared. Consequently, a variety of strategies have 
been adopted to discover Class sources: (1) Large scale imag- 
ing for outflows in the optical or near infrared (NIR; e.g. Ziener 
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& Eisloffel ll999l) : Stanke et al. 120021) ): ( 2) High resolution pro- 
cessing of IRAS data (e.g. Hurt & Barsony 1 1996); O'Linger et al. 
1 19991) ): (3) s ub-mil limetre or millime tre con tinuum mapping (e.g. 
Shirley et al. l2000l) : Motte & Andre <200ll) ) ; (4) d eep radio con- 
tinuum or line surveys (e.g. Bontemps et al. ll995T) : Bourke et al. 
Il997l) ). Therefore it is evident to keep in mind that source samples 
are subject to strong selection effects. 

We require a statistically significant observational sample of 
protostars with accurate individual properties. The latter demand a 
uniform observational coverage of the e ntire SED from the NIR to 
the millimetre range. Froebrich 120051) recently presented a cata- 
logue with all known confirmed and candidate Class sources in 
the literature. This catalogue contains 50 objects that possess suffi- 
cient observational data to compute the three main source proper- 
ties (bolometric temperature and luminosity, envelope mass; T bol , 
Lboi, M clw ) accurately. 

We can potentially test a range of cloud models. Each model 
assumes a particular environment, with initial conditions, bound- 
ary conditions, external for ces an d feedback effects to be consid- 
ered (Mac Low & Klessen (2004)). However, in the present study, 
our attention is focused on a set of accretion rates obtained from 
numeric al sim ulations of turbulent clouds (see § 12.21 Schmeja & 
Klessen 120041) : hereafter SK04). The mass distribution of the cores 
so formed can be compared to the initial mass function and the total 
mass fraction accumulated in protostellar cores ca n be related to the 
over all star- formation efficiency (Klessen 1200 lbl) . Padoan & Nord- 
lund ^20021) 1. These recent numerical studies are able to provide an 
almost complete model of star formation in clusters, inclu ding th e 
formation of disks and binary systems (see also Bate et al. (2003)). 
As a more fundamental test of these models, we need to compare 
model predictions to the observable properties of stars caught in the 
act of forming (Tboi, Lboi, M env ). The numerical models used here, 
as well as the works from other groups (e.g. Bate et al. 120031) ') are 
only able to follow the density and velocity of material. They are 
hence very useful in obtaining final mass spectra, binary fraction, 
etc.. However, they are not able to predict observables like bolo- 
metric temperature and/or luminosity. To derive those quantities, in 
principle a 3D radiative transfer calculation needs to be performed 
at each time step. Since this is much too time consuming we applied 
an evolutionary scheme to obtain those quantities (see Sect. l2.3> . 

Alternatively, mass accretion rates can be calculated from the- 
oretical and semi-empirical approaches. In the so-calle d "stan dard 
solution" of the collapse of an isothermal sphere (Shu ll977lf ) the 
mass accretion rate is constant at a value of 0.975 c] / G, where 
Cj is the sound speed and G the gravitational constant. Lar ger bu t 
also constant accreti on rates (47 c] / G) are found by Larson 1 19691) 
and Penston 1 19691) for the co llapse of initially uniform isother- 
mal spheres. Henriksen et al. 1 19971) find that mass ejection and 
mass accretion both decline significantly with time during proto- 
stellar evolution. Us ing a P lummer-type density profile Whitworth 
& Ward-Thompson 1200 ll) predict an accretion rate that possesses 
a very early peak during the Class phase and falls off rapi dly af- 
ter the object develops into a Class 1 protostar. Shu et al. 1 2004) 
obtained similar results through analytical and numerical works. 
Hydro dynamical s imula tions of star formation (Foster & Chevalier 
ll993l) : Tomisak a 1 19961) ) controlled by supersonic turbulence (e.g. 
Klessen 1200 id) : SK04) yield similar behaviour. 

In theory, there are two means to interface model and observed 
data sets. Ideal would be to derive a set of mass accretion rates 
directly from observations of protostellar cores but infall speeds 
are extremely difficult to measure. Hence, conversely, we convert 
the mass infall rates from the models into the observable quantities. 



The main three are the bolometric temperature and luminosity, Tboi 
and Lboi, and the envelope mass, M cnv . 

Several such conversion schemes have been develope d. Evo- 
lutionary models w ere fir st discussed by Bontemps et al. j^99g) 
and Saraceno et al. Il996l) . Along these lines, Andre et al. l2000l) 
took an exponentially decline in both envelope mass and accretion 
rate to predict L|-, n |-M tracks. An analytical scheme was presented 
by Myers et al. Il998l). They assumed a mass infall rate match- 
ing the Shu solution at early times and then displaying an expo- 
nential fall off with time. The core which supplies the mass also 
provides the obscuration, thus fixing evol utionary tracks on a Tb i- 
Lboi diagram. Based on this work, Smith ( Il998tl2000tl2002l) : S98 
hereafter) presented an evolutionary scheme but adopted alterna- 
tive analytical forms for the mass accretion (exponential increase 
and power law decrease) which resem ble the results obtained by 
Whitworth & Ward-Thompson 1200 ll) or SK04. This evolution- 
ary model, founded on mass transfer between the envelope, disc, 
protostar and jets, has provided successful interpretations for vari- 
ous sets of observational data (T bo i, L bo i, M env , ou tflow l uminosity) 
of Cla ss and Class 1 objects (e.g. D avis et al. Il998l) : Yu et al. 
l200d) : Stanke l2000l) : Froebrich et al. 120031) ). Here, we adopt the 
S98 scheme but replace the analytical accretion rates with those de- 
rived numerically. We note that this scheme does not include any 
inclination effects that become very important during the Class 1 
phase for T bo i and L b() i determinations. We thus focus our work on 
Class objects. 

Our approach in this paper is as follows. We combine dif- 
ferent sets of mass accretion rates obtained from gravoturbulent 
mo lecular cloud frag mentati on calculations (gt-models; Klessen et 
al. 120001) : Klessen 1200 id) : AppendixlAt. as described in detail 
in SK04, and the evolutionary scheme (e-model) from S98 (see 
also AppendixlBT For each set of mass accretion rates, we deter- 
mine evolutionary tracks in the Tboi-Lboi-M elw parameter space and 
compare them with the distribution of the observed Class sample. 
A 3D-Kolmogorov-Smirnov test (KS-test, see also AppendixO 
yields a probability for the two distributions to be drawn from the 
same basic population (called agreement hereafter). The agreement 
parameter will thus determine how well a model is able to explain 
the currently available set of observations. 

We remark that two separate models (gt-models and e-models) 
are needed to obtain the model protostellar quantities that can be 
compared to the observational data. Thus, we are only able to de- 
termine how well the combination of both can explain the obser- 
vations. We are not able to decide which model is responsible for 
potential disagreement with the observations. However, the varia- 
tion of free parameters in the models and a search for correlations 
of the parameter values with the agreement will help to define con- 
straints and provide suggestions for improved model combinations. 

In Sect.|2|we describe the observational data and our models. 
Results and the discussion are put forward in Sect.|3| Details of 
the gravoturbulent and evolutionary models and the Kolmogorov- 
Smirnov test are described in AppendixeslAltolCl 



2 OBSERVATIONS AND MODELS 
2.1 Observational data 

Sever al sam ples of Class source s have been pub lished (e.g. Chen 
et al. 1 19951): Andre et al. 120001) : Shirley et al. l2000l) : Motte & 
Andre 1200 ll) ). The latest study, which combines all previous ones 
and computes the object properties uniformly from all published 
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data, is presented by Froebrich 120051) . There, a list of 95 con- 
firmed or candidate objects was compiled. Fifty of these sources 
possess sufficient observational data to be classified as Class or 
Class 0/1 objects and allow us to determine the three source proper- 
ties Tboi, Lboi, and M clw accurately. To ensure a reliable comparison 
of models and observations, we apply the following restrictions to 
this sample: 

(1) We omit sources that have distances larger than 500pc. 
This reduces the bias towards higher mass objects. 

(2) All objects in Taurus are underluminous compared to the 
other sources, considering Tboi and M cnv . Note that this does not 
only mean a low luminosity for these objects, but rather a low lu- 
minosity combined with a high envelope mass. The three very low 
luminosity objects in our sample (see Fig.0 combine low luminos- 
ity (1 L Q ) with low envelope mass (0.2 M Q ) and might form (very) 
low mass stars. In total, 25 % of the Class objects in the sample 
are underluminous. We exclude them for the comparison between 
the models and observations (but see the discussion in Sect. l3.6> . 

(3) A histogram of the number of sources in a certain Tb r 
bin shows that the observational sample is very incomplete at high 
bolometric temperatures. Hence the comparison of observations 
and models is restricted to Tboi < 80 K. 

This selection procedure leaves us with a sample of 27 Class 
sources. This sample consists mostly of objects in Perseus, Orion 
and Serpens. 

2.2 Gravoturbulent models 

We performed numerical simulations of the fragmentation and col- 
lapse of turbulent, self-gravitating gas clouds and the resulting for- 
mation and evolution of protostars. We employed a c ode b ased on 
smoothed particle hydrodynamics (SPH; Monaghan ll992l) ) in or- 
der to resolve large density contrasts and to follow the evolution 
over a long tim escale . The code includes periodic bou ndary c ondi- 
tions (Klessen j 19971) 1 and sink particles (Bate et al. <1995l) 1 that 
replace high-density cores while keeping track of mass, linear and 
angular momentum. A detailed study of the 24 sets of protostel- 
lar mass accretion rates resulting from the gt -models is presented 
in SK04 (see also Jappsen & Klessen fe004t) ) and we adopt their 
gt-model notation. 

We now define when a model core is considered to be a Class 
object. Observationally, a protostar i s a Cla ss source according to 
the d efinitio n given in Andre et al. feOOfJ) . It is shown by Andre 
et al. i 19931) that this definition is roughly equivalent to the phys- 
ical state where the mass in the envelope exceeds the mass of the 
central protostar. Thus, a model protostar is considered a Class 
object when (a) the mass of the central core exceeds 10~ 2 M G (this 
corresponds to the ti me whe n the second hydrostatic core is formed 
in the centre (Larson fc003l) ) and distinguishes the Class sources 
from the pre-stellar core phase) and (b) the mass of the envelope is 
larger than the mass of the central star. This separates the Class 
from the Class 1 objects. 

Small number statistics will always be a major concern when 
comparing models and observations of Class objects. In order not 
to introduce further uncertainties from the model side, we restricted 
our analysis to those gt-models that possess more than 37 stars and 
have a numerical resolution of at least 2-10 5 particles. This leaves 
16 out of the 24 sets of accretion rates of SK04 for further analysis. 
In addition, accretion rates for individual protostars from the gt- 
models are smoothed by a viscous time scale of the accretion disc 
of ~ 10 4 yrs (see AppendixlAlfor details). 

A small fraction of model protostars were highly accelerated 



Table 1. Variable parameters from the e-model of S98, together with the 
range they were varied in and the original value. The last two columns pro- 
vide the overall range found to lead to the best agreement. 
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1.0 
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100 


30 


30 


80 
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0.5 
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1.75 


1.4 


3.2 


R m [AU] 
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100 


30 


35 


80 


M eff [%] 





50 


30 


36 


44 
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1.4 


2.0 


1.5 


1.55 


1.80 


k [cm 2 g~'] 


2 


6 


4 


3.0 


5.0 



(e.g. by ejection from a multiple system). Due to the adopted pe- 
riodic boundary conditions in our calculations, these objects cross 
the computational domain many times while continuing to accrete. 
In reality, however, these protostars would have quickly left the 
high-density gas of the star-forming region and would not be able 
to gain more mass. We therefore consider accretion to stop after the 
object has crossed the computational box more than ten times. 

The radius of a sink particle is fixed at 280 AU, the physics 
(e.g. exact accretion, radiation) inside this volume cannot be re- 
solved. Therefore, an evolutionary model describing the processes 
inside the sink particle is required. 



2.3 Evolutionary scheme 

The e-model is used to transform the mass accretion rates from the 
gt-models into observable quantities such as Tboi, L DO i, and M en y It 
is based on mass and energy transfer between the different compo- 
nents of the forming star (protostar, disc, envelope, jet). Mass con- 
servation is the main principle in this analytical model. Each simu- 
lated accretion rate from the gt-models is taken as the mass inflow 
rate from a spherical envelope onto the inner disc/protostar/jet sys- 
tem. All the mass flows through the disc but only a fraction accretes 
onto the protostar. The rest is ejected within two jets. The ejected 
mass fraction is small but not negligible in the Class stage. See 
AppendixlBlfor a detailed description of the parameters and equa- 
tions of the evolutionary model used in this work. 

Many parameters and constants are needed to fully describe 
the en velope and predict the radiative properties (Myers et al. 
Il998l) . S98, AppendixlBl. We carefully chose a subset of param- 
eters which we kept variable: (1) r env ; The temperature of the 
outer envelope, out to which the envelope mass is determined. (2) 
/raCen V ; The total mass that will be accreted onto the protostar is 
distributed in the envelope and the disc. A constant ratio of these 
two masses is assumed and frac em represents the fraction of the 
total mass that is in the envelope. (3) M cxtra ; This is the supplemen- 
tary mass fraction in the immediate surroundings which does not 
actually fall towards the star, to be accreted or jetted away, but is re- 
moved directly from the core probably through a feedback process. 
See Eq. lB7l for the exact definition used. (4/5) to, a; These two pa- 
rameters describe the dispersion of the extra mass in the envelope. 
(6) R m ; The inner radius of the envelope. (7) M c ff ; The maximum 
percentage of the infalling envelope mass which is ejected into the 
outflow. (8) p; The power-law index of the density distribution in 
the envelope, assumed constant in time. (9) k; The opacity of the 
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envelope material at 12yum. In TableQwe list the parameter ranges 
which were tested. 

We now only need to impose a mass accretion rate from the 
gt-models for the evolution of a model protostar to be fully deter- 
mined. Using all individual accretion rates in a gt-model allows us 
to build up a T bol -L bo i-M cnv distribution which can be compared to 
the observational dataset. 

Some of our free parameters can be, and have been, observed 
for individual objects. This includes for example the power-law in- 
dex of the de nsity distribution in the envelope (e.g. Chandler & 
Ric her feOOCl) . Motte & Andre <200ll) 1 or frac cm (e.g. Looney et 
al. 120031) 1. We chose to keep those parameters variable in order 
to test our method, since we expected to obtain values within the 
observational constraints. 



2.4 3D KS-Test and probabilities 

The best way to compare two distributions of data points is through 
a Kolmogorov-Smirnov test. This test yields the probability of two 
distributions being drawn from the same basic population. In one 
dimension, the KS-test compares the cumulative probability func- 
tions of a sample and a model and for large sample sizes the prob- 
ability can be determined analytically. Since our distributions are 
three dimensional (T bo i-L bo i-M clw ), there are no analytical means to 
perform such a test. We therefore gen eralise d the method of a two 
dimensional KS-test (e.g. Singh et al. 12004) ) for our purpose. The 
basic result of this test is the value D 3D that ranges from zero to 
one. The lower this value the better the two distributions match. 
Using a Monte Carlo method £> 3D can be converted to the agree- 
ment (P3D), which gives the probability that the two distributions 
are drawn from the same population. 

The task of finding the parameter set for the e-models that re- 
sults in the highest agreement is a multi-dimensional, non-linear, 
minimisation problem. We solved this using a Monte Carlo ap- 
proach. In Appendixlclwe outline the details of this method which 
allows to constrain the range for the e-model parameters and to de- 
termine the best agreement. 



3 ANALYSIS & DISCUSSION 
3.1 Evolutionary tracks 

One of our goals is to investigate the accuracy to which the final 
mass of a protostar (M flm i) can be estimated from its present loca- 
tion in the T bo i-L bo | diagram. To achieve this, we sort for each gt- 
model the individual model stars into final mass bins with a width 
of 0.3 in logarithmic units (equivalent to a factor two in mass). The 
following ranges for the mass bins were chosen: <0.2, 0.2... 0.4, 
0.4...0.8, 0.8... 1.6, 1.6.. .3.2, >3.2M Q . The particular size of the 
mass bins was adopted in order to ensure a reasonably large num- 
ber of stars in each bin for the gt-models. Depending on the mass 
function and the total number of objects in the gt-model, there are 
up to 20 objects per bin for the lower masses (M flmll < 1.6 M Q ). The 
higher mass bins naturally suffer from a paucity of objects. 

We determined for every gt-model the average accretion rate 
history in each mass bin. These mean accretion rates are then used 
to determine evolutionary tracks in the T bol -L bo ] diagram. The left 
column of Fig.Q shows the tracks for the six mass bins of the 
models M05k8, M2k2, and M6k2a as an example. The lowest fi- 
nal masses correspond to the thinnest line, and so forth. Note that 
the average evolutionary tracks show higher L bo i at a given T bo i 



for higher final star masses. This represents the fact that stars with 
higher final masses on average have higher accretion rates (SK04) 
and hence bolometric luminosities. 

Further we analysed the individual evolutionary tracks. As an 
example, we plot all individual tracks for the same three models 
with final stellar masses in the range 0.8 to 1.6 M G (mass bin 4, 
solar mass stars) in the middle column of Fig.Q These tracks are 
determined assuming typical e-model parameters from the range 
that best fit the observations (see below). We find that all tracks 
show a similar general behaviour. The tracks start off at T bo i « 15 K 
for solar mass stars. Lower mass Class sources start off at higher 
temperatures. Then they show a gradual increase in luminosity until 
the end of the Class phase. Some notable exceptions are evident 
in Fig.Q which imply that T bo i is not always a reliable guide as to 
the envelope-protostar mass ratio. 

Are we able to estimate the final mass of a Class source from 
its position in the T bo i-L bo i diagram? To address this question, we 
determined the 1 cr scatter of the individual tracks at each time step 
in all mass bins. In the right column of Fig.Qthis scatter is shown 
for mass bin 4 as thin lines. Almost independent of the gt-model, 
we find that the scatter has about the size of the separation between 
tracks for two adjacent mass bins. Hence, we conclude that from a 
certain position in the T bo i-L bo i diagram, we are only able to esti- 
mate the final mass of the protostar to within a factor of two. Note 
that this estimate does not take into account possible errors in the 
measurement of T bo i and L bo i, as well as possible different accretion 
histories. This supports the notion that other observables might be 
more adequate fo r datin g Class protostars, e.g. the L smm /L bo i ratio 
(Young & Evans 120051) ). in agreement with th e origi nal observa- 
tional definition of Class sources (Andre et al. Il993h ). 

We further find that tracks determined from averaged accre- 
tion rates noticeably differ from the averaged evolutionary tracks 
of the individual stars. This is evident in the right column of Fig.Q 
which shows as a thick line the track from the averaged accretion 
rate. At some points, it approaches or even lies outside the ± 1 a 
range of the individual evolutionary tracks. This clearly shows, to- 
gether with the findings in the above paragraph, that the individual 
accretion history significantly influences the position in the T bo] - 
L bo i diagram, in addition to the final stellar mass. 

The general behaviour of the T bo i-L b „i tracks (characteristics 
of the individual tracks, the average tracks, and the scatter around 
the average) is independent of the considered mass bin and gt- 
model. 

3.2 Distribution in T bol -L bo i-M cnv 

In Fig.|2|we show how the model tracks (M05k8, M2k2, M6k2a, 
from left to right) are distributed in the full T bo i-L bo i-M elw parame- 
ter space. The grey-scale background of the individual panels gives 
the probability to find one of the model stars at this particular posi- 
tion. Darker regions indicate higher probabilities. We used typical 
e-model parameters (see Table[2) to determine these diagrams. For 
comparison the observational datapoints are overplotted. Note that 
we plotted only model points with T bo i < 80 K, as this is the limit 
of the observational data and only these are used in the comparison 
with the observations. 

As evident in Fig.|5| the model distributions cover roughly the 
same area as the observations. Especially in the T bo i-L bol and T bor 
M cnv plane the models are able to explain the peak in the observed 
distribution. In the L bo i-M ellv plane, however, the models fail to ex- 
plain this peak (see right column in Fig.[3J. Our models predict 
smaller envelope masses (by a factor of two) compared to parts 




Figure 1. Evolutionary tracks in the Tboi-Lbol parameter space of the model Class sources in the gt-models M05k8 (top row), M2k2 (middle row), and 
M6k2a (bottom row). The tracks are determined using typical e-model parameters from the range best fitting the observations (see Tablel2l. Left column: 
Tracks for the average accretion rates in the six mass bins. Thicker lines correspond to higher final masses. Note that there are no stars in mass bin 6 in the 
M6k2a model. Middle column: All individual tracks in mass bin 4. Right column: Track for average accretion rate of stars in mass bin 4 (solar mass stars; thick 
line) and the ± 1 <x scatter of the individual tracks (thin lines). The vertical lines mark ages of the sources starting from 5 ■ 10 3 yrs in steps of 5 ■ 10 3 yrs. The 
open diamonds in each panel mark the positions of the observational sample (taken from Froebrich 2005). See Sect. l3.1l for more details. 



of the observations. Considering the cluster of observational points 
at about M cnv = 1 M , the discrepancy could be a selection effect 
in the observations. In addition, observed en velope masses are un- 
certain by a factor of two (Motte & Andre |200l|)). See also the 
discussion in Sec . 13. 71 



agreement with the estimated final stellar mass spectrum in the co r- 
responding observational sample (r = -0.9±0.2; Froebrich 120051) ). 
Recall, that for this sample also no binary correction has been 
made. 



3.3 Initial mass function 

We analyse the mass function of the final masses of model stars 
(IMF) in the gt-models in order to compare it to that of the obser- 
vational sample. The protostellar mass functions of the gt-models 
show a decline for masses larger than about 0.3 M Q . The mean value 
of the slope of all considered models is (r> = -0.84 in the mass 
range -0.5 < log M/M Q < 1. Th is is less steep than the Salpeter 
slope of -1.35 (Salpeter 119551) ') but our results are biased by the 
fact that binaries or multiple systems cannot be resolved. At the 
low-mass end, the mass function is constrained by the SPH reso- 
lution limit of ~ 0.05 M . Nevertheless, the mean value is in good 



We test if the slope of the IMF in the gt-models is related to 
the best agreement (P30). However, no systematic dependence of 
P3D on T is present. In particular, a good agreement of T between 
the gt-model and the observations does not necessarily lead to a 
high /> 3D -value. As discussed in the last paragraph (see also the 
middle panel of Fig.0, the differences in the position in the T bor 
Lboi diagram are partly due to the accretion history and not due 
to the final mass. Hence, the test performed here is much more 
sensitive to the distribution of individual accretion histories than to 
the slope of the IMF. 
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(see Tablel2l. Only model points with T bo i smaller than 80 K are shown in the plots since 



3.4 Evolutionary model 

To constrain the free parameters in the e-model, we investigate here 
how the chosen values influence the agreement P 3D . Table|2|con- 
tains the range of values which provide the highest agreement for 
each gt-model. In particular we show the values for the best fitting 
5 % of the e-models. It is more meaningful to provide these ranges 
instead of the particular parameters for the single best fitting model. 
These are statistically less significant, since not the whole parame- 
ter space is tested by our method. We obtain the following results 
for all gt-models: 

• The temperature at the outer envelope boundary is lower than 



the standard value of 24 K. Typically temperatures between 15 and 
19 K lead to a good agreement. 

• A wide range for inner envelope radii (R m ) leads to high P^. 
Values from 35 to 80 AU are typical. This wide range is in agree- 
ment with the fact that a change in R in , while fixing all other pa- 
rameters, usually results in very small changes of />3d- However, 
the range for the best fitting models does not include the hitherto 
standard value of 30 AU. 

• In most cases small amounts of extra mass in the envelope 
(1.0 < M cxlra < 2.4) lead to high agreement. 

• The parameter a shows a wide range of possible values for a 
good agreement (1.4 < a < 3.2). A similar wide range is found for 



Evolution of Class protostars: Models vs. Observations 1 



Table 2. Parameter ranges needed in the e-model to best match the observational data (see AppendixIBlfor details on the parameters). We list the values for 
the best fitting 5 % of the parameter combinations. The last two columns list the median and average duration of the Class phase of the individual stars in the 
models. 

Model T em frac em M exlra to a R m M eff p k D 3d P 3d ^ t^ Q 

[K] [%] [10 3 yrs] [AU] [%] [cmV'] W [10 3 yrs] [10 3 yrs] 



G2 


15 


19 


81-89 


0.6-2.4 


35-75 


1.4 


3.2 


30 


80 


44-47 


1.6-1 


9 


3.0-5.2 


0.30-0.37 


37.2- 


5.93 


41 


44 


M01k2 


14 


18 


87-95 


1.0-2.8 


20-85 


1.4 


3.4 


40 


75 


28-44 


1.6-1 


8 


3.0-4.6 


0.37-0.41 


12.2- 


2.43 


39 


41 


M05k4 


17 


20 


89-94 


1.2-2.0 


35-75 


1.6 


3.2 


45 


85 


41-44 


1.5-1 


7 


3.0-4.4 


0.29-0.37 


40.3 - 


4.72 


31 


36 


M05k8 


16 


19 


90-95 


0.9-2.0 


25-80 


1.9 


3.0 


40 


85 


40-46 


1.5-1 


7 


3.0-4.6 


0.36-0.43 


45.9- 


6.53 


31 


34 


M2k2 


15 


20 


86-92 


1.0-2.4 


35-85 


1.6 


3.2 


35 


80 


35-43 


1.5-1 


7 


3.4-5.2 


0.29-0.36 


73.1 - 


19.6 


32 


39 


M2k4 


15 


18 


86-92 


1.0-2.4 


35-75 


1.4 


3.2 


40 


75 


39-43 


1.5-1 


7 


3.2-5.2 


0.48-0.52 


30.5 - 


4.69 


15 


23 


M2k8 


13 


18 


80-90 


1.0-2.2 


30-75 


1.2 


3.6 


35 


75 


32-41 


1.5-1 


7 


3.4-5.4 


0.27-0.34 


67.4- 


17.0 


45 


56 


M3k2 


15 


18 


82-94 


0.6-2.6 


20-90 


1.4 


3.6 


35 


85 


37-44 


1.5-1 


8 


2.8-4.6 


0.28-0.35 


55.1 - 


10.8 


44 


43 


M3k4 


12 


14 


80-92 


0.6-1.6 


25-80 


1.4 


3.4 


30 


70 


32-44 


1.5-1 


S 


3.2-5.2 


0.33-0.40 


5.99- 


0.53 


87 


105 


M6k2a 


16 


20 


85-95 


1.0-2.4 


20-85 


1.4 


3.4 


45 


85 


35-43 


1.6-1 


8 


2.8-4.6 


0.32-0.36 


46.5 - 


20.4 


26 


32 


M6k4a 


15 


20 


93-97 


0.8-3.2 


40-80 


1.4 


3.2 


30 


70 


12-34 


1.6-1 


8 


3.2-5.2 


0.47-0.50 


0.05 - 


0.01 


38 


120 


M6k2b 


15 


20 


82-94 


1.2-3.0 


35-80 


1.6 


3.2 


35 


75 


36-40 


1.6-1 


8 


3.2-5.4 


0.35-0.40 


14.2- 


2.84 


26 


35 


M6k2c 


15 


20 


86-96 


1.0-3.0 


35-75 


1.4 


3.2 


35 


75 


36-42 


1.6-1 


9 


3.4-4.8 


0.38-0.43 


5.30- 


0.62 


29 


51 


M6k4c 


13 


18 


82-85 


0.6-2.6 


20-85 


1.2 


3.2 


30 


80 


28-44 


1.5-1 


8 


3.0-5.2 


0.45-0.51 


22.2- 


2.57 


14 


54 


M10k2 


16 


22 


92-96 


1.0-3.6 


25-70 


1.2 


3.0 


35 


80 


10-42 


1.6-1 


9 


3.2-5.4 


0.60-0.62 


0.003- 


0.00 


19 


24 


M10k8 


14 


17 


93-98 


0.8-3.2 


30-85 


0.8 


3.8 


35 


80 


16-38 


1.5-1 


7 


3.4-5.2 


0.43-0.49 


2.60- 


0.23 


44 


124 



to which ranges from 30 to 80 • 10 3 yrs. The larger values found for 
t a suggest that the additional mass stays longer in the envelope to 
sustain low bolometric temperatures. 

• A range of values of p is found to be able to best explain the 
data. However, in the majority of the cases 1.55 < p < 1.80 leads 
to the best results. These values are in between the theoretical so- 
lutions for a free falling envelope (1.5) and a singular isothermal 
sphere (2.0). The range found here seems to be a compromise be- 
tween younger and older objects. Since there are probably more 
older sources in the sample, the best values are closer to p = 1.5. In 
principle an evolution from p = 2.0 to p = 1.5 would be expected. 
Our findings do not contradict this. Our obtained values are in good 
agreem ent wi th (sub)-mm obse rvatio ns of envelopes (Chandler & 
Richer feOOOl) : Motte & Andre J200lh t. 

• We found the mass fraction in the envelope to be in the range 
from 86 to 96 %. This would imply a disc mass of 4 to 14 % of 
the envelope mass, mostly smaller than the standard value of 13 % 
used so far, but in agreeme nt with millimetre interferometric obser- 
vations (e.g. Looney et al. 120031) ). 

• The amount of mass ejected into the jets is found to range 
from 36 to 44 %. This is a narrow range and somewhat contradic- 
tory to the observations and jet launching models which predict 
that at maximum 30 % of the material is ejected into the jets. We 
interpret the high amount of ejected material here as a combination 
of two things. (1) Material really ejected into the jets. (2) Material 
accreted onto the star, but the resulting luminosity escapes directly 
through the cavity generated by the jets and is not observed. If the 
accretion luminosity is radiated uniformly away, then the additional 
ejected material would imply opening angles between 30 and 60 de- 
grees for the cavities, not an untypi cal val ue for many of the older 
Class sources (e.g. Padman et al. 1 19971) "). To conclusively prove 
this hypothesis one needs to include the outflow luminosity for all 
objects into the comparison of models and observations. 

• The best-fitting opacity range (3.0 to 5.0cm 2 g~') agrees very 
well with the standard value of 4cm 2 g~'. 



Q 




10 30 50 70 90 110 

Mean Class duration [10 3 yrs] 



130 



Figure 3. Mean duration of the Class phase of the model stars vs. the 
P3D value for the gravoturbulent models. Circles mark models that possess a 
median duration of the Class phase between 2 and 6 ■ 10 4 yrs and a ratio of 
median to mean duration of the Class phase of higher than 0.7. Triangles 
mark models were one of these conditions is not fulfilled. Gravoturbulent 
models with Mach numbers in the range 0.5 s? M < 3.0 are shown as filled 
symbols, models outside this range as open symbols. The Gaussian collapse 
model is marked as 'G'. The dashed lines enclose the region where a high 
agreement of the models with the observations was found. 



3.5 Gravoturbulent models 

We now investigate if the initial conditions of the gt-models in- 
fluence the agreement with the observations and the best fitting e- 
model parameters. The first point to note is that we find no signif- 
icant correlations of any of the obtained parameter values of the 
e-models with the Mach number or wavelength that characterise 
the turbulence in the gt-model. 

There are, however, large differences in the quality of agree- 
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ment between the various gt-models and the observations. In the 
following we investigate the cause of these differences. We deter- 
mined for each set of accretion rates the median and the mean du- 
ration of the Class phase of the model stars. The median for all 
gt-models ranges from 14 to 87 • 10 3 yrs, while the mean duration 
spans from 23 to 124 • 10 3 yrs. The much wider range for the aver- 
age duration is due to some model stars (outliers) that go through 
an extremely long Class phase. These objects have a large, sec- 
ondary accretion peak which in many cases might not be physical 
but is caused by the lack of feedback mechanisms. Due to the mass 
criterion used, such peaks shift the transition from Class to 1 to 
later times. In Table|5|we list the median and mean duration of the 
Class phase of all stars in the models in the last two columns. 

All gt-models with P 3D -values larger than 30 % are associated 
with a mean duration of the Class phase between 2 and 6 • 10 4 yrs. 
Most of them possess a ratio of median to mean duration of the 
Class phase exceeding 0.7. This indicates that no or only a few 
outliers are present in the set of accretion rates. When applying 
these two conditions to all gt-models, we select nine of them (cir- 
cles + G in FigQ. Only M01k2 (12.2 %) and M6k2b (14.2 %) do 
not possess P 3 d > 30%. Hence 80% of the gt-models with these 
properties lead to high agreement. There are seven models that do 
not fulfil one of the criteria (triangles in Fig.[3). Six of these possess 
P 3 D-values smaller than 30 %. Hence 85 % of the gt-models that do 
not comply with both criteria lead to low agreement. This shows 
that high _P 3D -values are only obtained for gt-models that generate 
few outliers and the 'right' duration of the Class phase. 

The duration of the Class stage shows a linear correlation 
with the median of the fit parameter t of SK04. This parameter in- 
dicates the time it takes until the peak accretion rate is reached. It is 
related to the local free-fall time and, thus, can be used as a rough 
estimate of the average density of the core at the onset of collapse 
(see SK04 for details). Given our best fitting models and the corre- 
sponding Class lifetimes, we infer a range of 1-5 10 5 cnT 3 for the 
local density at the start of the collapse. This is in good agreement 
with measured central densities of prc-stellar cores (e.g. Bacmann 
et al. l2002h ). Our results thus suggest that stars in the observational 
sample may only form in pre-stellar cloud cores that reach a central 
density of order of 10 5 cirr 3 . Obviously, the duration of the Class 
phase depends directly on the initial local density. 

As can be seen in Table|2| gt-models with Mach numbers 
0.5 < M < 3 best describe the observations (filled symbols in 
Fig-0}- Six out of seven such models possess a high P 3D -value, 
while only one out of eight models with M > 6 has a high 
P 3D . The exceptions M3k4 (M < 6, P 3D =5.99%) and M6k2a 
(M > 6, P 3 d=46.5 %), as well as the fact that similar gt-models 
(e.g. M6k2a/b/c) result in very different P 3 D-values are, however, 
understandable since protostellar collapse is, in essence, a rather lo- 
calised process, connected to the turbulent cloud environment only 
via the available mass and angular momentum inflow rate. It indi- 
cates that star formation is a highly stochastic process, influenced 
more strongly by local properties than by the global initial condi- 
tions, a fact which is also found, for example, by Bate et al. 1200 3l) . 

The range of Mach numbers found as best explanation for 
the data reflects the composition of the observational sample. It 
mainl y cont ains sources from Perseus, Orion, and Serpens (Froe- 
brich ( 2005)). There typical Mach numbers are in the range 1 < 
M < 6 (e.g. Castets & Langer ll995h . Schmeja et al. <2005h ). It 
is, however, interesting to note that the majority of gt-models with 
M > 6 are not able to explain the observations, independent of 
the choice for the e-model parameters. This indicates that the ini- 



tial conditions do influence the observational properties of Class 
protostars. 

3.6 Underluminous sources 

We recall that we excluded 25 % of the observational data in our 
determination of P 3 d- Those excluded exhibit a much lower bolo- 
metric luminosity than the other objects, for their given bolometric 
temperature and envelope mass. They could be objects which expe- 
rience either (a) quiescent states of low accretion or (b) a different 
time dependence of the accretion rate (Froebrich 120051) ). 

A detailed investigation of the individual accretion rates from 
the gt-models rules out the existence of quiet accretion phases as 
the dominant cause for the low luminosity of these sources. On 
average, the individual model stars spend less than 5 % of their time 
in a phase of suppressed accretion which we define as being a factor 
of six lower than the average value. This implies that we should 
observe only about 5 % of the sources in such a state, compared to 
25 %. Hence, these objects appear to follow a different accretion 
history. 

We investigate this in more detail by performing the procedure 
of determining P 3 d using not only the 'normal' sources, but also 
the Taurus-like objects. We find that this reveals worse agreement 
with P 3D values typically halved. The agreement between models 
and observations is always lower when we try to explain the Taurus 
like sources. This strengthens the proposal that in the case of 'nor- 
mal' sources, turbulence governs the accretion process, while other 
processes are responsible for the underluminous objects. Consid- 
ering their lower luminosity and hence lower accretion rates, we 
speculate that ambipolar diffusion might be an important process 
governing accretion in these sources. 

A more detailed analysis of a sub-sample of Taurus like ob- 
jects would be desirable since the time dependence of the mass 
accret ion rate might vary from region to region (Henriksen et al. 
Il997t) ). Even tho ugh m ore and more such objects are discovered 
(e.g. Young et al. l2004h ). indicating that the known percentage of 
these objects is clearly a lower limit, the current sample still suffers 
from too small number statistics to perform a reliable statistical 
analysis. 

3.7 Further discussion 

As evident from Table|2| the best agreement between observations 
and models is rather low (» 70 %). Here we will discuss three pos- 
sible causes for this: 

1. Our procedure: There are many free parameters which, in 
principle, could take values within a large range, as noted in Ta- 
bleQ A straightforward test of all possible combinations to solve 
this highly non-linear minimisation problem of many variables is 
impossible. We thus applied a Monte Carlo method to obtain the 
best fitting parameter combination (see Appendixlcl. The obtained 
probability P^> is a lower limit since not all possible parameter 
combinations have been tested. Note that only slightly smaller val- 
ues for Z) 3D can lead to probabilities much closer to 100 %, espe- 
cially for models that already possess a high agreement. 

When converting D 3D to P 3D we restricted the random selection 
to model values with bolometric temperatures smaller than 80 K. 
This is the limit of our observational sample of Class sources. 
There are further restrictions in the observational data. No objects 
outside a given range in L bo i and M clw . This is most likely an ob- 
servational bias. We refrain from applying the other observational 
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limits for the random selection process for the P 3D determination 
because the observational limits are not well defined/understood. 
If doing so, however, we would increase the probabilities for the 
agreement in many cases significantly towards 100 %. 

2. Observational data: As stated above, the observations 
clearly cover only a limited range of the T bo i-L bor M env parameter 
space. All limits, except T bo i < 80 K, are observational biases. Note 
that the observational sample consists of sources from several star 
forming regions. It is possible that e.g. the sensitivity limit for the 
bolometric luminosity is not the same in all regions, changing the 
observed statistics of low luminosity/mass sources compared to the 
real one. 

As shown in Fig.Q applying a fixed limit for T bo i as a divide 
between Class and Class 1 objects is not valid (see also Young 
& Evans 120051) ). Most of the solar mass stars undergo this tran- 
sition at temperatures close to 80 K. However, depending on the 
individual accretion history of each star, some quite significant 
deviations from this value are evident. Furthermore, lower mass 
Class 0/Class 1 transition objects tend to be warmer than the tem- 
perature limit applied here. Thus, the object sample might suffer 
from misclassified objects. This also shows that there is a differ- 
ence in the definition for Class protostars based on observations 
and models. Furthermore, the observational definition describes a 
different physical state of the object depending on its accretion his- 
tory and final mass. We overcome most of these uncertainties, how- 
ever, by applying the T bo i < 80 K restriction when determining />3d- 

3. The models: All our models are constructed on the most ba- 
sic principles (e.g. the gt-models neglect feedback mechanisms, the 
disc mass is a constant fraction of the envelope mass, etc.). In ad- 
dition, all model stars are treated with the same set of e-model pa- 
rameters. It is not clear if some of the parameters might depend on 
the final mass of the star. Many parameters are also kept constant 
in time. This might not be valid (e.g. the inner envelope radius, the 
fraction of mass in the disc, and the powerlaw index of the enve- 
lope density distribution could depend on the evolutionary state, 
final mass of the star). It is in many cases, however, uncertain how 
these dependencies can be parameterised. We thus kept these pa- 
rameters constant so as not to introduce even more free arbitrary 
parameters. 

Given these three points, we conclude that a very good agree- 
ment between the models and observations should not be expected. 
It would rather be very worrying if our method would result in a 
99 % agreement between models and observations. We further have 
to consider the possibility that star formation, and in particular the 
mass accretion process, might take such a diversity of paths that no 
simple unifying model can expect to capture them all, leading to 
low agreements even for future much improved observational sam- 
ples and models. 



4 CONCLUSIONS 

We extracted mass accretion rates within cores generated by gravo- 
turbulent simulations and inserted them into a simple protostellar 
evolutionary scheme to calculate protostellar evolutionary tracks. A 
principle dependence of the position of these tracks in the T bo |-L bo i 
diagram on the final mass is found. A detailed analysis, however, 
shows that we are not able to determine the final mass of a particu- 
lar object from its measured bolometric temperature and luminos- 
ity more accurately than a factor of two. The particular shape of 
an evolutionary track is largely determined by the accretion history 
and by its final mass. Hence, in the context of the gravoturbulent 



scheme, unique evolutionary tracks do not exist in the Class or 
Class 0/1 phase, as opposed to the ensuing pre-main sequence evo- 
lution. It also implies that T bol is not always a reliable guide to the 
envelo pe-protostar mass ratio, as suggested also by Young & Evans 
12005). 

A 3D-KS-test was used to compare the distribution of our evo- 
lutionary tracks in the T bo i-L bor M cnv parameter space with an ob- 
servational sample of Class objects. By varying free parameters 
associated with the evolutionary scheme, we are able to determine 
constraints for some of the parameters (r env , frac cm , M ex(ra , M cff , 
p). Others were found to have no significant influence on the agree- 
ment between models and observations (f , a, R m , k). A comparison 
of the parameter values obtained by our method with observational 
constraints from individual objects (e.g. for frac cm or p) can be 
conducted. This supports the reliability of our method as similar 
parameter values are obtained. 

Only gravoturbulent models generating model stars with a 
specific density at the start of the collapse and Class lifetimes 
between 2 and 6 • 10 4 yrs lead to a good agreement with the obser- 
vations. This Class lifetime is consistent with estimates of a few 
10 4 yrs based on num ber rati os or dynamical timescales of outflows 
(see e.g. Andre at al. <2000h ). It is shown that the Class lifetime 
is correlated with the local density at the onset of the collapse. The 
determined density range agrees well with density measurements 
for pre-stellar cores. 

Gravoturbulent models with 0.5 < Ai < 3 lead to the best 
agreement with the observations. The failure of models with M > 6 
to explain the data indicates that the initial conditions influence the 
observational properties of Class sources. The initial mass func- 
tion appears not to influence the agreement systematically. The 
highest probability that our distributions of model tracks and ob- 
servational data points are drawn from the same basic population is 
70 %. This is reasonable given the uncertainties in the observations 
and the simple assumptions in the models. However, the method 
can be readily adapted to compare larger future source samples, 
different sets of accretion rates or other envelope models. 

Applying the KS-test to both the 'normal' and underluminous, 
'Taurus-like' objects in the source sample, we find that the deter- 
mined probability P3D is about twice as large when testing the 'nor- 
mal' objects. This is consistent with the proposal that turbulence 
governs the accretion rates in the majority of the objects. 
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APPENDIX A: GRAVOTURBULENT MODELS 

Our simulations describe the evolution of (1) two globally unstable 
model clouds that contract from Gaussian initial conditions without 
turbulence and (2) 22 models where turbulence is maintained with 
constant root mean square Mach numbers. The models are labelled 
G1/G2 for the Gaussian, and MMkk (with rms Mach number At 
and wavenumber k) for the turbulent models, following SK04. The 
designators a, b, c distinguish models that have the same M and 
k values, but different random realisations of the turbulent driving 
field; they also differ in the time when self-gravity is 'switched on'. 

The gt-models are computed in normalised units. To scale 
to physical units, we adopt a mean density of n(H 2 ) = 10 5 cur 3 
and a temperature of 11.3K corresponding to a sound speed c s = 
0.2kms~'. For the two Gaussian models, the total mass present is 
220 M and the size of the cube is 0.34 pc. For the turbulent models, 
the cube contains a mass of 120 M Q within a volume of (0.28 pc) 3 . 
The mean thermal Jeans mass in all models is (Mj) = 1 M G and the 
global free-fall timescale is T ff = 10 5 yr. Requiring that the local 
Jeans m ass is a lways resolved by at least 100 gas particles (Bate & 
Burkert ll997l) 1. the resolution limit is 0.058 M Q in all our turbulent 
models and 0.044 M Q in the Gaussian models. 

The radius of a sink particle is fixed at 280 AU. Infalling gas 
particles undergo several tests to check if they remain bound to the 
sink particle before they are considered accreted. We cannot re- 
solve the evolution in the interior of the control volume defined by 
the sink particle. Because of angular momentum conservation most 
of the infalling matter will accumulate in a protostellar disc within 
which it is transported inwar ds by viscous and p ossibl y gravita- 
tional torques ( Prin gle ^98 lh : Papaloizou & Lin 1 19951) : see also 
Jappsen & Klessen 120041) for the models discussed here). The lat- 
ter will be provided by spiral density waves that develop when the 
disc becomes too massive. This occurs when mass is loaded onto 
the disc faster than it is removed by viscous transport. Altogether, 
the disc will not prevent or delay material from accreting onto the 
protostar for long. It acts as a buffer and smoothes possible accre- 
tion spikes. For the mass range considered here feedback eff ects ar e 
too weak to halt or dela y accre tion (Wuchterl & Klessen 1200 lh : 
Wuchterl & Tscharnuter 120031) ). With typical disc sizes of the or- 
der of a few hundred AU, the control volume encloses both pro- 
tostar and disc. The measured core accretion rates are hence good 
estimates of the actual stellar accretion rates. Strong deviations may 
be expected only if the protostellar core fragments into a binary sys- 
tem, in which case the infalling mass is distributed over two stars. 
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In addition, if material with very high angular momentum is ac- 
creted, a certain mass fraction may end up in a circumbinary disc 
and not accrete onto the star at all. 

The timescale over which accretion is smoothed by the disc 
can be approximated by the viscous timescale 

t y * a- y [ (R/H) 2 ts, (Al) 

(Pringle il98ll> t. with the viscosity parameter a v , the disc radius 
R, the disc scale height H, and the rotational period = Us- 
ing typical values of or v = 0.01, R/H = 10 and % = 1 yr this 
yields a rough estimate of t v ss 10 4 yrs for the viscous timescale 
at a radius of 1 AU. A large fraction of mass will lie external to 
this radius, leading to a larger f v , but this will be compensated for 
by stronger gravitational torques, which can lead to significantly 
larger values of the effective viscosity in the disc (see e.g. Laugh- 
lin & Rozyczka Jl99fih ). Since we cannot model the detailed be- 
haviour inside the disc, we take 10 4 yrs as a rough measure for the 
smoothing timescale. This value is used to smooth the individual 
mass accretion rates from the gt-models. Our results do not sig- 
nificantly depend on the smoothing scale that we apply. Adopting 
viscous timescales of 5 or 15- 10 3 yrs yields very similar results. 

There are further smoothing effects and time-shifts, e.g. the 
luminosity generated in the centre typically escapes a Class en- 
velope in about 300 yrs (as determined by a random walk model). 
Hence, this time lag is neglected in comparison to the disc vis- 
cosity timescale. The viscosity of the disc has additionally the ef- 
fect that material transported from the envelope onto the disc needs 
some time before it gets accreted onto the star and the accretion 
luminosity is generated. Hence, the observed envelope masses and 
luminosities do not represent the same time in the evolution. The 
observed luminosity corresponds to that of an envelope mass which 
was present about a viscosity timescale earlier. However, since this 
effect does not become significant until the disc is formed, it can be 
neglected in the very early stages of the evolution. In later stages, 
the evolution along the T bol -L bo i track is slower and the luminosity 
only slightly changes within one viscosity timescale. In particu- 
lar the change is smaller than typical errors of the measurements. 
Hence, we also neglect this effect in our calculations. 



APPENDIX B: EVOLUTIONARY MODEL 

Accretion rates applied in the e-model to date assume specific func- 
tional forms with a fixed time scale for all protostars. As a conse- 
quence, even without a statistical comparison, it fails to explain the 
data. To ove rcome this, new parameters were introduced by Myers 
et al. Il998l) . They took different evolutionary timescales for the 
envelope and the accretion, and subsequently found that the obser- 
vations required (1) the accretion time scale to exceed the envelope 
timescale by a factor of a few and (2) the final mass of the star to 
be quite a small fraction of the original envelope mass. 

In addition, there are several free parameters in the e-model. 
Some of their values are not as yet tightly constrained through 
observations or theory. Plausible values, however, predict evolu- 
tionary tracks across the same r egion as occupied by the observed 
protostars (e.g. Froebrich et al. 120031) ). Nevertheless, an immedi- 
ate problem encountered with the basic model was that the proto- 
stars spent little time in the early Class phase but lingered in the 
late Class phase for periods exceeding 10 5 yr, implying contrary 
statistics to those observed. To overcome this, a distinct mass com- 
ponent was introduced into the envelope. In the e-model, this arbi- 
trary mass would not accrete, but would be dispersed directly. The 



rationale was that star formation remains inefficient due to outflow 
activity from the protostar. The second arbitrary parameter was the 
fraction of mass contained within a flattened distribution. That is, it 
is assumed from the outset, that a small fraction of the gas is avail- 
able to accrete but does not add to the obscuration of the protostar. 
However, given these additions, it becomes difficult to relate the 
predictions uniquely to the chosen parameter set through the KS- 
test if we attempt to calculate analytically which parameter values 
would lead to the highest agreement between models and observa- 
tions. 

Additional free parameters in the e-model are an inner density, 
and inner and outer radius for the envelope. We have generalised the 
previous studies by introducing a power-law index (p) for the radial 
density distribution, where logp oc -pxlogR, rather than fixing the 
value at 1.5. Instead of the inner density we employ the total enve- 
lope mass as the variable. In addition, rather than the outer radius, 
we adopt the outer temperature r cnv as the variable which then de- 
termines the outer radius. Most of the mass initially contained in 
the envelope will eventually fall onto the protostar. The accretion 
luminosity depends on the collapse radius and the contraction of the 
central hydrostatic object. Here, we assume that a constant density 
core develops. 

Equations 

In the following we describe in more detail the equations and pa- 
rameters used in the evolutionary model for this work. 

The initial mass of the envelope and circumstellar disk, M , 
is found by integrating the mass accretion rate, MJf), over time. 
Hence, the envelope/disk mass is given by 



M e ,(f) = M( 



o-J M,(f 



)-df. 



The mass just in the spherical envelope is taken as 
M cm (t) = frac em ■ M c , d (r). 



(Bl) 



(B2) 



In the first approximation, the accreting material is taken as instan- 
taneously reaching the central hydrostatic core, defining the proto- 
stellar radius, R,. A fraction (1 - rj{t)) of this material is accreted 
onto the star i.e. 



M,(t) 



t 



(1 - r 1 {t'))-M^t')-dH. 



(B3) 



The parameter rj(t) will depend on the jet launching mechanism. 
We suppose that it is a function of the relative accretion rate i.e. 



j](t) = M ei 



M a (f) 



A/n, 



(B4) 



where M max is the peak accretion rate and £ = 2 has been assumed 
to date. Thus, M c ff is the maximum fraction of material ejected into 
the jets. 

We assume here that the hydrostatic core grows at constant 
density during the very early stages of star formation (no energy 
released through contraction): 



R« 



4.24 ■ 



M,(t) 


1/3 


M " 


2/3 


I M e 









(B5) 



We have found that alternative treatments produce similar results. 
The particular dependence on envelope mass corresponds to the 
case where the final escape speed of the jets is constant. 
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The gravitational energy released is assumed to end up in ra- 
diation or in the jets. We thus write the accretion luminosity as 



= (1 - n(t)) ■ G ■ MM ■ M,(t)IR,(t). 



(B6) 



The initial cloud core mass could exceed the envelope mass. 
That is, some extra mass is introduced. This mass is dispersed 
rather than being accreted. In other words, we can check to see if 
there is some form of feedback or turbulent dissipation in operation 
which maintains a cooler core during the early stages. We assume 
a form 

\-2<n 



M corc (f) = Me nv (f) 



1 +M„, 



t + tg 

to 



(B7) 



where a and t describe the dissipation of M exlra . Note that for 
A^extra = the particular values for to and a are superfluous and 
M core (t) is equal to M cnv (f). 

The radius corresponding to the outer edge of the envelope is 
determined by an assumed outer temperature, r e nv: 



1 



2;r ■ cr 



I/: 



(B8) 



where cr is the Stefan-Boltzmann constant. 

Given an inner radius for the envelope, (possibly fixed by 
the centrifugal barrier condition) and a power-law radial distribu- 
tion of density with index -p, the inner density of the envelope can 
be written 

(3 - p) ■ M COK (t) 



Pin(0 = 



An-Rl-e{t) 



where 



8{f) 



( Xw(0 
R m 



3-J> 



1. 



The inner temperature of the envelope is 

1/4 



Zta(f) 



In-cr-R- 



(B9) 



(BIO) 



(Bll) 



This yields an optical depth through the envelope proportional to 
the emissivity: 



■PmO)-Rk 



1 - e(ty-p 



1 



(B12) 



Following Myers et al. |199q). we take k as the emissivity at 12yum, 
A = 1.59 1(T 13 cm 2 g _1 Hz" 1 , h as the Planck constant, k the Boltz- 
mann constant, and calculate the bolometric temperature: 



Zt»i(r) 



r(9/2) ■ f (9/2) 

r(5) • f(5) 



h ■ k ■ T m (t) 



k-A-r e (t) 



1/2 



(B13) 



APPENDIX C: 3D KS-TEST AND PROBABILITIES 

To evaluate how well observations and models match, we need to 
determine the KS-parameter D^- First, we note that each data point 
(T! ii ., MJ, nv ) can be used to dissect the parameter space into 
eight octants. The difference D\ between the fraction of model data 
points and observational data points in each of the eight octants 
in the T bo i-L bor M env parameter space is calculated for each data 
point. Here, the index ( covers both the observational and model 
data points. In case of N data points, this procedure results in &N 
numbers D\. The maximum absolute difference value is then de- 
fined as the KS-parameter D 3D = MAXflAI). Note that this value 
can range from zero to one. 



As described in Sect. l2.ll and l2.2l the restrictions in the source 
sample and models lead to the problem that we have to com- 
pare two samples with extremely different sizes. The observational 
sample consists of 27 points. Depending on the set of accretion 
rates we use, the model sample consists of about ten thousand 
points (100 tracks of individual stars and on average 100 points 
per single track, considering our timestep of 300 yrs and assum- 
ing a lifetime of 3-10 4 yrs for Class sources). Hence, the test 
described above has to be performed about 10 4 times to deter- 
mine Z?3d for one case, resulting in a huge amount of computa- 
tions. Since we are only interested in the maximum of the abso- 
lute differences between the relative numbers of the model points 
and observational points in the octants, we can constrain the in- 
vestigated area to a cuboid containing all the observational dat- 
apoints. A test outside these boundaries will not lead to abso- 
lute differences larger than obtained inside the cuboid. In other 
words, MAX(|Z);|) corresponds to a point situated inside or at 
the boundary of the cuboid defined as: [T^' mm < T bo , < T* s ' max ; 

Lobs.min „ t , t obs.max -« *obs,min , » if ^ a *obs,maxn tt - 
bol < L bo i < L bol ; M env ' < M cnv < M env ' ] . However, since 

most of the model points are still situated within this cuboid, the 

required amount of computations is not significantly reduced. 

Hence, we do not perform the test at all data points within this 
area, but rather at a certain number of points inside the cuboid. We 
performed extensive tests in order to establish the best compromise 
between computing time (number of points where the test has to be 
performed) and the accuracy of the method. As a result, we found 
that a 5 x 5 x 5 point grid in the cuboid is sufficient to determine £>3d 
with the same accuracy as the probability of the agreement between 
observations and models (see below). The ranges in T bo i, L bo ], and 
M env are thus divided into four equally sized regions. While for 
T bo i linear spacing was chosen, L bo i and M cnv are divided into equal 
spaces in logarithmic units. 

In order to convert the £>3d value into a probability of agree- 
ment P 3D , a Monte Carlo method has to be applied. For this pur- 
pose, we generated artificial observational data points by randomly 
selecting 27 model points out of all model points and treating them 
as observational data. Since the observational data is limited to 
T bo i < 80 K, we applied this limit for the random selections as 
well. Then the KS-test was performed (using the same grid posi- 
tions as with the real observational data points), leading to a value 
D^§- Repeating this process, a distribution N(£>^-) is created. The 
probability can then be determined by: 



J N(< c ) 



3 D 



1 

/ 



N(< c ) 



(CI) 



All £>3o-values need to be converted to P3D by such a Monte 
Carlo simulation. To minimise the computational needs we tested if 
this is needed for every single £> 3D -value. It turned out that the dis- 
tributions N(Z)^-) obtained by a Monte Carlo simulation for one gt- 
model but different parameter sets in the e-model are similar. Hence 
all Z) 3D -values of one gt-model can be converted into P 3D using the 
same distribution. However, distributions N(£>^) for distinct gt- 
models are partly very different. Thus, to find the best fitting model 
combination we selected the parameter combination leading to the 
smallest £>3D-value for each gt-model to determine N(D^J-). This 
was then used to calculate the agreement P 3D . 

The total number of random selections to build up the distri- 
bution N(Z)^") was a critical point concerning the accuracy of the 
method. Hence, we tested for which number of repeats the final 
values of P iD did not vary by more than one percent if the same 
test was repeated. For our distribution of model points, this was the 
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case for 3-10 4 random selection processes. This accuracy is not de- 
graded by the 5x5x5 grid, chosen for performing the KS-test, 
instead of doing the test at all datapoints in the cuboid. 

There are nine free parameters in the e-models which we kept 
variable. Our task here is to find the combination of these param- 
eters that leads to the best agreement with the models. A test of 
all possible combinations with an accurate spacing requires an im- 
mense amount of computations. We thus decided to use a Monte 
Carlo approach for this task as well. We randomly selected 100 
values for each parameter from the range given in TableQand de- 
termined £>3d for these 100 random combinations. We then selected 
the parameter sets that lead to £> 3D -values which were smaller than 
(Z)"™ + Z)"™)/2. The maximum and minimum for each parame- 
ter value used in these selected sets was then taken to re-set the 
ranges out of which parameter values are selected. Then again 100 
e-models with randomly selected parameters out of the new range 
were tested. This procedure was repeated 30 times, finally leading 
to 3100 tested e-models. These models are used to determine the 
best fitting parameter ranges (Table|2j for the individual gt-models. 
These ranges were then used to determine another 10000 random 
parameter combinations, in order to find the best agreement. 

We tested for a subset of parameters and smaller ranges if the 
Monte Carlo method leads to the same results as a straight forward 
testing of all possible parameter combinations. Both, the smallest 
Z?3D-value and the best fitting ranges for the parameter values, are 
found to be similar. Note that the smallest found D3D-value has to 
be considered an upper limit, since not the whole parameter space 
was tested. Alternatively the probability is a lower limit. 



